//################################################################################
// MIT License
// Copyright (c) 2024 ZhangYihua
//
// Change Logs:
// Date           Author       Notes
// 2023-03-20     ZhangYihua   first version
//
// Description  : 
//################################################################################

module fft_sim #(
parameter           IDW                     = 12,       // data width
parameter           IFW                     = 9,        // fractional width
parameter           N                       = 8192,     // must be 2^N_BW
parameter   string  WIN                     = "NONE",   // "HANNING"
parameter           PLOT_EN                 = 1'b1,     // plot spectrum on or off
parameter   real    CUT_OFF                 = 0.125,    // ideal cut-off frequency
parameter   real    ROLL_OFF                = 0.25,     // filter roll-off factor, range [0,1]
parameter           PRINT_EN                = 1'b1      // print infor on or off
) ( 
input                                       rst_n,
input                                       clk,
input                                       cke,

input       signed  [IDW-1:0]               dat_i,      // s(IDW, IFW)
input       signed  [IDW-1:0]               dat_q,      // s(IDW, IFW)

input                                       act,        // each toggling H/L start a transform
output  reg                                 done
);

//################################################################################
// define local varialbe and localparam
//################################################################################
localparam          N_BW                    = $clog2(N);
localparam          PI                      = 3.1415926536;
localparam          TAB_NUM                 = N/2;

real                                        di_ary[N-1:0];
real                                        dq_ary[N-1:0];
real                                        ci_ary[TAB_NUM-1:0];
real                                        cq_ary[TAB_NUM-1:0];
real                                        win_ary[N-1:0];
real                                        pwr_ary[N-1:0];

real                                        fft;
integer                                     freq;

real left_max, right_max, mid_avg;

//################################################################################
// main
//################################################################################

initial begin:TAB
    integer         i;
    real            ph;
    
    for (i=0; i<TAB_NUM; i=i+1) begin
        ph = (-2.0)*PI*i/N;     // WN(k) = e^(-j*2*Pi*k/N)

        ci_ary[i] = $cos(ph);
        cq_ary[i] = $sin(ph);
    end

    if (WIN=="HANNING") begin
        for (i=0; i<N; i=i+1) begin
            ph = PI*(2*i+1)/N;
            win_ary[i] = 0.5 - 0.5*$cos(ph);
        end
    end
end

initial begin:SMP_FFT
    integer i, j;
    real phs, ri, rq, pwr;

    fft <=`U_DLY 0;
    freq <=`U_DLY 0;
    done <=`U_DLY 1'b0;

    @(posedge rst_n);

    forever begin
        @(posedge act or negedge act);

        // store and resort I/Q samples
        for (i=0; i<N; i=i+1) begin
            @(posedge clk);
            while(cke==1'b0) begin
                @(posedge clk);
            end

            ri = dat_i/(2.0**(IFW));
            rq = dat_q/(2.0**(IFW));
            j = invert_order_f(i);
            if (WIN=="NONE") begin
                di_ary[j] = ri;
                dq_ary[j] = rq;
            end else begin
                di_ary[j] = ri*win_ary[i];
                dq_ary[j] = rq*win_ary[i];
            end
        end

        // butterfly round
        for (i=1; i<=N_BW; i=i+1) begin
             butterfly_t(i);
        end

        // to dB
        for (i=0; i<N; i=i+1) begin
            ri = di_ary[i];
            rq = dq_ary[i];
            pwr = ri*ri + rq*rq;
            if (pwr==0) begin
                pwr = 1.0/(2.0**(IFW+16))/(2.0**(IFW+16));
            end
            j = (i<N/2) ? (i + N/2) : (i - N/2);
            pwr_ary[j] = 10*$log10(pwr);
        end

        analyse_t;

        if (PLOT_EN==1'b1) begin
            plot_t;
        end

        done <=`U_DLY 1'b1;
        @(posedge clk);
        done <=`U_DLY 1'b0;
    end
end

function [N_BW-1:0] invert_order_f;
    input [N_BW-1:0] seq;

    integer i;
    begin
        for (i=0; i<N_BW; i=i+1) begin
             invert_order_f[i] = seq[N_BW-1-i];
        end
    end
endfunction

task butterfly_t; 
    input integer round;

    integer grp, pair;
    integer i,j,s;
    real ai,aq,bi,bq;         
    real ci,cq;
    real mi,mq;
begin
    grp = 1<<(N_BW-round);
    pair = 1<<(round-1);
    for (i=0; i<grp; i=i+1) begin
        s = i*pair*2;
        for (j=0; j<pair; j=j+1) begin
            ai = di_ary[s+j];
            aq = dq_ary[s+j];

            bi = di_ary[s+j+pair];
            bq = dq_ary[s+j+pair];

            ci = ci_ary[j*grp];
            cq = cq_ary[j*grp];

            mi = bi*ci - bq*cq;
            mq = bi*cq + bq*ci;

            di_ary[s+j]      = ai + mi;
            dq_ary[s+j]      = aq + mq;

            di_ary[s+j+pair] = ai - mi;
            dq_ary[s+j+pair] = aq - mq;
        end
    end
end
endtask

task plot_t; 
    integer i;

    for (i=0; i<N; i=i+1) begin
        @(posedge clk);
        fft <=`U_DLY pwr_ary[i];
        freq <=`U_DLY i - N/2;
    end
    @(posedge clk);
    fft <=`U_DLY 0;
    freq <=`U_DLY 0;
endtask

task analyse_t;
    integer i, k, m;
    integer fcut, froll;

    fcut = CUT_OFF*(N/2);
    froll = ROLL_OFF*fcut;
    i=N/2-fcut+froll;
    k=N/2+fcut-froll;
    mid_avg = 0;
    m = 0;
    while(i<=k) begin
        mid_avg = mid_avg + pwr_ary[i];
        m = m + 1;
        i=i+1;
    end
    if (m>0) begin
        mid_avg = mid_avg/m;
    end

    left_max=pwr_ary[0];
    i=1;
    k=N/2-fcut-froll;
    while(i<=k) begin
        if (pwr_ary[i]>left_max) begin
            left_max = pwr_ary[i];
        end
        i=i+1;
    end

    right_max=pwr_ary[N-1];
    i=N-2;
    k=N/2+fcut+froll;
    while(i>=k) begin
        if (pwr_ary[i]>right_max) begin
            right_max = pwr_ary[i];
        end
        i=i-1;
    end

    if (PRINT_EN==1'b1) begin
        $display("%m Left Max:%.1fdB; Mid Avg:%.1fdB; Right Max:%.1fdB", left_max, mid_avg, right_max);
    end
endtask

//################################################################################
// ASSERTION
//################################################################################

`ifdef CBB_ASSERT_ON
// synopsys translate_off


// synopsys translate_on
`endif

endmodule
